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Abstract. A two-step optimization is proposed to represent an arbitrary quantum 
state to a desired accuracy with the least number of gaussians in phase space. The 
Husimi distribution of the quantum state provides the information to determine the 
modulus of the weight for the gaussians. Then, the phase information contained in 
the Wigner distribution is used to obtain the full complex weights by considering the 
relative phases for pairs of gaussians, the chords. The method is exemplified with 
several excited states n of the harmonic and the Morse oscillators. A semiclassical 
interpretation of the number of gaussians needed as a function of the quantum number 
n is given. The representation can also be used to characterize Wigner and Husimi 
distributions directly which do not originate in a quantum state. 
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1. Introduction 

In many situations there is a need to represent a given quantum state or a probability 
distribution of states in terms of a simple basis which can be easily handled in further 
processing, e.g., for integration in matrix elements. Gaussians form such a basis 
although it is overcomplete. Gaussians have the additional advantage that they can 
adopt the form of coherent states in phase space, thereby supplying a direct link to 
classical mechanics. Hence, the formulation in terms of gaussians is also very useful 
for a mixed quantum-classical or quasiclassical description of a system. In particular 
the evolution of systems with many degrees of freedom can often be described within 
classical molecular dynamics only. Here, one is interested in obtaining optimal initial 
conditions by translating the initial quantum state of the system into a suitable classical 
phase space distribution which can be classically propagated p. While all quantum 
information on an amplitude level is lost in this approximation, we will introduce a 
representation which can preserve basic quantum interferences by construction even in 
the case of the purely classical propagation. Finally, it is of general and basic interest 
to reconstruct amplitudes or the density matrix, i.e., coherent information, from phase 
space distributions. This has been done recently in terms of the moments in time 
of the Wigner distribution |2]. Although there is no principal need to do so we will 
concentrate on the Wigner distribution due to the attention it has received recently in 
different areas. In several experiments the Wigner distribution of microscopic systems 
has been measured jS]. Beyond the quantum context, in the time-frequency phase space 
of signal processing, analysis of the Wigner distribution j3] is an important tool (4j, as 
for optics speech analysis [H] or the monitoring of machinery conditions [7]. 

The decomposition of a given state |\E') into coherent states can be regarded as a 
particular form of phase space representation. Labeling positions by g or Q and momenta 
by p or P, we translate the ground state of the harmonic oscillator by X = (P, Q) to 
obtain the coherent state in position representation, i.e. its wavefunction is [HI E] 

For simplicity, we chose unit frequency (cu = 1) and mass (m = 1) for the harmonic 
oscillator without loss of generality. (Often the complex variable Z = {Q + iP)/\/h is 
used and the constant phase —PQ/2% is omitted). Even though the coherent state basis 
is overcomplete, the exact representation, 

m = - I dX\X){X\^), (2) 

is unique |H1 E]. However, such a decomposition into coherent states is monstrously 
inefficient, as shown by the overlap of coherent states. 




^^^v[-—^ + ^x^x■\, (3) 
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so that 

l(^|X')r = expf-^^^V (4) 



We have used the wedge product, 

X AX' = PQ' - QP' = (JX) • X', (5) 

to shorten the notation. (With the second equation we also define the symplectic matrix 
J). It follows that one can approximate any given state by a finite expansion 

|vI/)-^a,|X,) (6) 

j 

with arbitrary accuracy. One possibility would be to place the coherent states on a grid 
and then determine the complex coefficients aj by minimizing the error, defined in some 
convenient way. Of course, it is much better to use some features of |\E') as a preliminary 
guide to where the \Xj) should be located. 

An obvious necessary condition for a good fit is that the coherent state intensity 
|(X|\I')p in (jH) be well approximated by 

\{Xm'=\j:aAX\X,)\' 
j 

= E|a,n(X|X,)r + 5:aX(X|X,)(X|X,)* (7) 

for each phase space point X. But note that, if the diagonal terms in this expression 
are only appreciable for Xj very close to X, a nondiagonal term decays as a product of 
gaussians centered on a pair of separate phase space vectors. Such a term is negligible 
for all X unless Xj is very close to X^. Moreover, such nondiagonal terms are not 
positive definite, which justifies our performing the preliminary fit 

p^(X) = -\{Xm' ^ E \a,\'-\{X\X,)\' = Y: G,(X), (8) 



where 



G,(X)=exp(-^^^). (9) 



The factor l/vr has here been introduced so that the left hand side becomes the Husimi 
function 0, i.e., the Husimi representation of the state l^'), also known as the Q-function 
in quantum optics jTTj. Here, is then decomposed into gaussians phase space packets 
(gpps) which are Husimi representations of the coherent states, specified by 

The Husimi function in its own right uniquely defines a quantum state, though 
part of the quantum phase information is locked up into its delicate analytic properties. 
These are hard to translate into other representations, but we propose to sidestep this 
problem. Our method is to make a preliminary fit of the Husimi function (jH)), which is 
real and positive definite, and therefore much simpler than fitting the complex coherent 
states directly in (jH)). However we can now identify the coefficients of the Husimi fitting 
as being proportional to the square modulus of the coefficients of the superposition of 



Optimal representations of quantum states by gaussians in phase space 



4 



coherent states The phases of the coefficients are then determined by relating the 
Husimi function to the coarse-graining of the Wigner function, i.e. we bring in yet a 
third phase space representation of the state It should be pointed out that, once the 
state has been fitted by a finite number of coherent states, it can be viewed in any 
representation. The special attention given here to the Wigner function is justified by the 
way that it highlights quantum interferences. In fact such interferences are sometimes 
mimicked by extra gaussians in the context of numerical work with classical densities pp. 
Thus, in the next section we present the Wigner function, its relation to the Husimi 
function and its specific form for an arbitrary superposition of gaussians. In addition, we 
discuss and establish the differences between the representation of quantum states such 
as the wavefunction, the Wigner distribution and the Husimi distribution using these 
gaussians. The Wigner distribution sets the stage for the delicate procedure of fitting 
the phases, i.e., to determine the appropriate interference fringes between the gaussians 
as detailed in section 3. In this section we propose a method for determining the 
phase of the expansion coefficients by Fourier analysing the Wigner function. Actually, 
this resulting chord function or characteristic function (in quantum optics JUl) can 
be derived directly from the initial wave function, without obtaining first the Wigner 
function. The Husimi function, on the other hand, supplies a robust overall picture of 
the quantum states which can be used to place the gaussians. This will be discussed 
in section 4 using explicit examples for illustration. Emphasis will be put on general 
symmetry criteria for placing the gaussians and numerical demonstrations based on 
Monte Carlo fitting of sampled wavefunctions. Finally in section 5, we semiclassically 
interpret the results of section 4, namely the increase of gaussians needed to represent 
quantum states ipn of a system with increasing quantum number n. A summary in 
section 6 concludes the paper. 

For simplicity, all our formulae are presented for systems with a single degree of 
freedom, but they are easily generalized to high dimensional systems. 



2. Wigner functions and their relations 

The state is uniquely described by the real function of phase space points x = {p, q), 

wAx) = J ^(9 + ^l^)(^l9-^)e^p(-^M') (10) 

known as the Wigner function ^2]. Inversion of this symmetrized Fourier transform 
leads back to the wavefunction, while its intensity results from the projection 

\{qm' = J dpW^ix). (11) 

Though the Wigner function can be negative for some points x, its projection is always 
positive. Actually, this is just a particular case of the general relation, 

= 27rn f dxW^{x)W<i,{x), (12) 
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for arbitrary states 1$) and \^). Thus, introducing the Wigner function for the coherent 
state 

= _Lexp(-(i^). (13) 

we obtain the Husimi representation of \^) as 

= 2h I dxWx{x)W^{x) 

71 J 

^^^) W^ix). (14) 

This shows that the positive definite Husimi function results from the gaussian 
smoothing of the Wigner function. Inserting the coherent state expansion given by 
(jni) into the Wigner transform (fTn|) . we obtain the well known form of the corresponding 
Wigner function, 

W^ix)=Y.\a,\^Wx,+2Re \J2 (^A^jkix)] , (15) 

where we define the crossed Wigner function (or Moyal bracket jl3j ) 

W,u{x) ^ J ^(9 + ^l^.)(^;^k-^)exp(-^M') (16) 

{x -Xjkf\ ( if .1 




exp -- - X A 5Xjk + -Xk A Xj 



Here we use the abbreviated relations Xjk = (Xk + Xj) /2 and 6Xjk = Xj—Xk, as well as 
the skew product Comparing (fT3j) with ((Tj) we find the same gaussian form for the 
diagonal terms of the Wigner function as for the Husimi function, but the offdiagonal 
terms are no longer negligible. For the Wigner function, any pair of coherent states 
determines a third gaussian centered halfway between their two centres, of comparable 
amplitude to the diagonal contributions of the individual gaussians. This interference 
term is only cancelled by the smoothing employed to obtain the Husimi function, because 
of the phase oscillations. Their wavelengths in phase space are 27ih/\6X\, so that 
the Husimi function is increasingly damped with growing separation of these gpps. 
This is an excellent illustration of the way that the Husimi function highlights the 
classical structure in a quantum state, while it hides quantum phase information. This 
question will be further developed in section 4. Of course the phase information is still 
contained in the small offdiagonal terms in (|7j), because the Husimi function is fully 
quantum mechanical, but this information is certainly more accessible in the Wigner 
representation. 

We can now understand the differences between using the position representation, 
i.e. the wave function, the Wigner representation or the Husimi representation, to fit a 
given state by a set of coherent states. The position representation has the advantage 
that there are no interference terms between pairs of coherent states. However, if the 
state is highly oscillatory, such as a highly excited eigenstate of a Hamiltonian, this 
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needs to be fitted by coherent states with large momenta P, which also have narrow 
oscillations. Clearly the fitting procedure becomes very unstable with respect to small 
errors in the initial function if large momenta are necessary. In contrast, the diagonal 
terms in the expansion of the Wigner function are smooth gaussians distributed over 
phase space, but again we obtain increasingly narrow oscillations for the interferences 
around the midpoint of each pair. Of course, the oscillations of the wave function and 
the Wigner function are not unrelated, since the latter have the same wavelength as the 
wave intensity in the case of a pair of gaussians located at (±-P, Q). 

Finally, the Husimi function washes out oscillations, the more effectively the tighter 
they happen to be. We shall review in section 4 how the smooth phase space distribution 
singles out the basic classical structure within the quantum state, which can be stably 
fitted by real gpps. 

The cost is that only the location of the gaussians in phase space and the square 
modulus of their coefficients are then determined. In the next section we present a 
method for retrieving the phases of the coefficients by taking the Fourier transform of 
the Wigner function. 

3. Fitting the phases 

Both, the direction and the wavelength of the fringes that modulate the interference 
peak halfway between two gpps at Xj and Xk in the Wigner representation, are entirely 
determined by SXj^. However, these fringes are translated by changing the phases of 
the coefficients aj and a^, which we wish to determine. Clearly, the best way to recover 
these phases is to take the Fourier transform of the Wigner function: 

It is also possible to express this directly from the wave function 

W^iO= /rfg(g + ||vl/)(^|g-|)exp(-^epg), (18) 

which resembles the Wigner transform (^UJ. This function is known as the Woodward 
ambiguity function in communications theory ^1] or the characteristic function in 
quantum optics [TU]. Usually the sign of is reversed in the definition of W, but 
here we follow where it is simply called the chord function. A justification for this 
term appears immediately if we consider the case of a superposition of gaussians that is 
our concern : 

where 

*'.,«) = ^exp(44.Y,A{) (20) 
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^jkii) = —7=— exp 



jk) 



\2' 



and 



X exp (^(e - 5X,fc) A X,, + j j • ^^^^ 

Thus all the diagonal terms collapse onto gaussians centered on ^ = 0, whereas each 
interference term is a gaussian centered on the chord joining Xj and X^. (Both directions 
±5Xjfc are present, leading to symmetric contributions in chord space. Indeed, we must 
have W{—^) = W*{^)). Let us now suppose that a given chord 6Xjk is sufficiently far 
from all the other chords, so that at ^ = SXj^ the chord function is dominated by the 
single nondiagonal contribution W{6Xjk), i.e. 

W^{6X^k) - E \ai\'Wu{6X^k) ~ a,a*iy^.,(5X,,) 
I 

aja* i Xk A Xj 



exp(-^- — ^). (22) 



2V27r^ h 
Then, defining 

aj = \aj\exp{i9j) (23) 

and recalling that we have already fitted \aj\ from the Husimi function, we obtain 
6jk = 6j — 6k from 

e,k - ^X, A X, = phase{W^{6Xjk) - E \ai\'Wu{6X,k)). (24) 



Of course, the modulus of both sides of (|22|) should be approximately equal, which 
provides a check on the previous Husimi fitting. The diagonal sum that is substracted 
from Wif,{6Xjk) in decays exponentially with SXjj^, so it will only affect the 
phases of the smaller chords. It might appear to be only consistent with our previous 
approximation of the Husimi function to neglect the diagonal sum in (j221); but this 
would perturb the phases of small chords. 



4. Placing the gaussians: examples 

So far, nothing has been said about how to place the gaussians in phase space that 
are meant to approximate a given Husimi function. An obvious criterion is that any 
relative maximum of this smooth distribution should also receive a coherent state. 
Further knowledge of the state should also be used. For instance, if it is known to 
be the eigenstate of a given Hamiltonian H with energy En, then the semiclassical 
considerations in the following section allow the restriction of the gaussian centers to 
\H{X) — En\ < h. Of course, the fitting procedure is only useful for cases in which H is 
unknown, or else |\E') may itself be a linear superposition of states, such as result from a 
two slit experiment. Here we immediately feel the advantage that the decomposition (jH]) 
is defined in terms of the state itself, rather than the density operator |\I/)(\E^|. Thus, the 
linearity of |\E') = ail^E'i) + 021^2) allows us to fit |\E'i) and |\E'2) independently and then 
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to superpose the fitted state as = ai\ipi)+a2\4'2)- A further simphcation results from 
possible symmetries. For instance, for the eigenstate of an even potential, one should 
place the gaussians symmetrically at {Pj,±Qj) with aj- = ±aj+, so as to guarantee 
even or odd states. Again, this is a very particular case, but time reversal invariance 
is much more common. If the state |\E') resulted from initial real wave functions and 
the evolution proceeded through interactions that preserve time reversal symmetries (in 
usual practice, if there are no magnetic fields) then the final state can also be real. In 
this case, one must choose symmetric pairs of gaussians at {±Pj,Qj) with equal real 
coefficients, to ensure that |\E') is also real. This halves the number of 

independent coherent states to be fitted and reduces the choice of phase to either zero, 
or 71. Any further knowledge should generally be used to reduce the randomness of the 
positions of the coherent states to be fitted to |\E'). The semiclassical considerations in 
the following section may be a further guideline. When this knowledge is exhausted, the 
best course is to optimize random guesses by Monte Carlo fitting of arbitrary gaussians, 
within the given constraints. 

Specifically, in order to measure the quality of the fitting, we define the relative 
error 



as the mean square deviation between the known Husimi function p^{x) of a quantum 
state \ip) and its fitting to the superposition of gpps defined in Q. This error a has to 
be minimized for a set of M grid points xi = {pi,qi)- To fit this superposition of gpps 
as closely as possible to Ptp{x), we proceed as follows: (i) the centers Xk = {pk, (lk)k=i,..,N 
of the coherent states are obtained by Monte-Carlo sampling the region where the 
distribution is significant; that is the region where the distribution exceeds a certain 
threshold S; (ii) the centers must not be too close, this means a minimum distance A 
between centers is fixed a priori, thereby avoiding overlaps and reducing the total number 
of gpps used for the representation; (iii) starting with = 1 we increase N by one 
and repeat the process until the desired accuracy is reached. Since only the coefficients 
Ak are fitted, minimizing the error cr reduces to a problem of linear optimization which 
is equivalent to solving the real matrix equation 



where the phase space points xi = {pi,qi)i=i,..,N that define the matrix elements 
Gik = Gk{xi) = exp — -^^^^f^) are randomly selected. Here, A = {Ai, A2^ A^) 
contains the coefficients to be determined and = (pi,P2, ■■■,Pn) those of the known 
Husimi function. In general, the matrix G is expected to be sparse and a variety of 
special algorithms can be used to solve (j^UI) efficiently j2I|. In principle one may also 
consider to optimize this fitting with respect to the widths of the gpps; but since the 
Husimi function is always smooth and positive definite, it is convenient to keep them 
frozen. In what follows, we have tested our method on the fitting of the eigenstates 
of the harmonic oscillator (HO) and that of the Morse oscillator (MO) supporting 18 




(25) 



GA = p^ 



(26) 
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Figure 1. Husimi distributions of the 3rd excited state of the harmonic oscillator (HO) 
and for the 5th excited state of the Morse oscillator (MO); (top : original, bottom: 
fit). The fitting of HO (left) is reproduced with N = 7 {S = 0.035, A = 0.5) and that 
of MO (right) with N =U{S = 0.02, A = 1.3). 

bound states \TT\. We have not used the knowledge that these states have time reversal 
symmetry in the way mentioned earlier, so we should obtain results of comparable 
quality even for states with complex wave functions. Figure 1 shows, on the top, the 
known Husimi function for the 3rd excited state of HO (left) and for the 5th excited 
state of MO (right), with parameters u = m = h = 1. On the bottom the corresponding 
fitted distributions are shown. We found that = 7 is sufficient to reproduce quite 
well the Husimi function of HO whereas 14 coherent states (gpps) are needed in the 
MO example, to achieve a global relative error of a < 0.01. The respective values 
for the threshold 6 and the minimum distance A are given in the caption of figure 1. 
This preliminary fit on the Husimi function provides us with the total number of 
gpps, their location in phase space {pk, qu) and the expansion coefficients A^. These 
coefficients are related to those of the coherent states and their phases 9k by 

ak = y^27rMfcexp(i6'fc). (27) 

Up to this point, the phases 6k are still unknown. We only need to determine the 
phases differences 9jk in (fTK|) which can be easily extracted from (j211). This says that, 
for gpps, A^(A^ — l)/2 chords {9jk) are required to reproduce the Wigner function. 
These N{N — l)/2 pairs of gaussians (chords) are symmetrically spaced about ^ = 0, 
so that, being more dense, they overlap more than the original set of gaussians in 
phase space. Figure 2 shows on the top the Wigner function of the 3rd excited state of 
the HO together with its corresponding fit obtained with 7 gpps. The present fitting 
reproduces the Wigner function everywhere in phase space, even for negative regions 
which are signatures of quantum interference. This is clearly illustrated in figure 3 where 
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5 S 

Figure 2. Wigner distributions of the 3rd excited state of HO (top) and its 
corresponding fit (bottom) with = 7 as in figure 1. 




Figure 3. Deviations of the Monte Carlo Wigner sampled of the 3rd excited state of 
HO (left) and of the 5th excited state of MO (right). 

the deviations of the Monte Carlo sampled Wigner function, less than 10~^ everywhere 
in phase space, are depicted for both, the HO (left, n = 3) and the MO (right, n = 5). 
The global relative error a is subsequently found not to exceed 0.01. Moreover, we 
have computed the number N of gpps needed to represent the Wigner functions of the 
eigenstates of our test systems, as can be seen in figure 4. One sees that N grows 
proportionally to ^/n as the vibrational quantum number n increases. In the following 
section, we provide a semiclassical justification of this result by means of the WKB 
quantization. 

5. Semiclassical approach 

A semiclassical state 1^^)30 is supported by a curve in phase space, in the case of one 
degree of freedom, or in general by a surface with half the phase space dimension [THl llTj. 
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2 4 6 8 10 12 14 16 18 

vibrational quantum number n 

Figure 4. Number of coherent states N used for the fitting versus the vibrational 
quantum number n for both, harmonic oscillator (<)•, dashed line: fit with N = 
2.77V^+ 2.10) and Morse oscillators (• , solid line: fit with N = 4.87^^+2.51). 



If the curve or surface is closed, it must be Bohr-quantized and it will be symmetric about 
p = in the case of time reversal symmetry. In the case of the position representation, 
for each branch of the curve (surface), Pj{q), one defines the action 



The amplitudes aj{q) can also be expressed in terms of the actions, but the important 
point is that they are purely classical so that the only /i-dependence occurs in the 
exponential. Thus locally, for any small range of positions, the semiclassical wave 
function reduces to a superposition of plane waves characterized by the wave-vectors 
Pj{q)/h. The semiclassical Wigner function is also defined in terms of an action Sj{x) 
with respect to the classical curve (surface). But instead of the area between the curve 
and the g-axis, we are now concerned with the area sandwiched between the curve and 
one of its chords. The latter, C,k{x), is selected by the property that it is centred on the 
point X, as shown in figure 5. Thus, for one degree of freedom jlHI, 




(28) 



leading to the generalized WKB wavefunction. 




(29) 




(30) 
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Figure 5. Geometrical illustration of the chord £,k{x) in phase space x{p,q). 



with straightforward generahzations ^H]- The important point is that ^3] 
dSk 



dx 



(31) 



where we use the symplectic matrix J defined in (0). Thus the semiclassical state is 
again represented by a superposition of waves, 

'5xA^fc(x) 7r\ 
1 iJ ' 



W^,{x)sc ~ Yl COS ( - 

k \ 



(32) 



but now these are phase space waves. Comparison of this expression with (|22j) 
immediately reveals that these waves have exactly the same wave vector as a pair 
of gaussians at the tip of each chord ^fc(x). It follows that an important feature of 
a semiclassical state is automatically reproduced by fitting it with gaussians precisely 
placed on the corresponding classical curve (surface). Notice that, in fact, this is a major 
a priori obstacle to performing the fitting, because the oscillations near the midpoint 
of a long chord are very fine for % « Sk- One would need very narrow gaussians 
indeed to fit these phase space waves directly, rather than having them arise naturally 
as interferences. All that is left to determine the phase: S^/h — 7r/4. How can we 
single out the region in which to place the gpps, if we have no a priori knowledge of the 
classical structures correponding to a quantum state? The obvious course is to smooth 
the Wigner function with a Gaussian window. In other words, we initially fit the Husimi 
function. Recalling that the gpps have linear width of order a/^, while the wavelength of 
the Wigner oscillations are h/\^k\^ it follows that this gaussian window erases effectively 



Optimal representations of quantum states by gaussians in phase space 13 

a semiclassical Wigner function, except in the limit of very small chords. This is a simple 
explanation of the well known fact that only in the neighbourhood of the classical curve 
(surface) itself is the semiclassical Husimi function appreciable. In the simple case of 
one degree of freedom, the Husimi function is concentrated near the energy shell. Thus, 
by fitting gpps to a Husimi function, one automatically samples its relevant classical 
mainifold, if it happens to have a (possibly unknown) semiclassical structure. A basic 
assumption in the above argument is that the gpps along the curve (surface) are neither 
too crowded, to avoid confusing superpositions, nor too sparsely spaced, lest gaps should 
arise in the fitting. It is thus easy to estimate the growth of the number of coherent 
states required to fit the eigenstate \n) with the degree of excitation. Since the linear 
width of the gaussian scales as a/^, the phase space area of the curve grows as nh and 
hence its length grows as \fnE, the number of equispaced gaussians needed to cover 
the curve grows as y/n. Of course, we have assumed here that the shape of the curve 
does not change with n, as for the harmonic oscillator. If the eigencurve elongates for 
higher excitation, the number of gpps necessary for a good fit will grow as ra", with 
1/2 < a < 1. We can also estimate the growth in time of the number of gpps required 
to fit an evolving semiclassical state with arbitrary precision. According to the theory 
of van Vleck [12], it is sufficient to evolve the curve (surface) classically and then to 
reconstruct the wave function. This principle can also be applied to Wigner functions, 
according to Berry and Balazs [201 ■ If the driving Hamiltonian is chaotic, then the 
curve will stretch at a rate depending on the Lyapunov exponent, A. Clearly, this also 
determines the initial rate of growth of the number of gpps needed for an adequate fit. 
Ultimately, when the curve covers densely all the available phase space (the energy shell 
of the driving Hamiltonian), a steady state will be reached where the number of gpps 
saturates. 

6. Conclusions 

Optimal fitting of quantum states by phase space gaussians is achieved by first fitting 
the Husimi function and then determining the quantum phases from the chord function, 
i.e. the Fourier transform of the Wigner function. It is clear from the discussion in 
the introduction that the density of gaussians in phase space must be finely adjusted: 
If ther are too few, essential features of the state will be missed, whereas an excessive 
number of gaussians would introduce interference terms in the Husimi function itself, 
which could only be accomodated by a much more complex variation of this method. 
The numerical examples in section 4 indicate clearly that , as well as achieving excellent 
overall accuracy, all essential qualitative features of the states are captured by this 
method, using a basis of gaussians that grows more slowly than the excitation number 
of the fitted states. This last result and further insight into the fitting procedure, follows 
from the semiclassical analysis in the previous section. This can be generalized readily 
to states of quantum systems with higher degrees of freedom if these are eigenstates 
of integrable systems (i. e. if they are supported by a Lagrangean surface in phase 
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space; see e. g. |I7j). Of course, more gaussians will then be needed for the fitting, 
but the preliminary fit of the Husimi function should still provide optimal results. The 
potential for this method to deal with the eigenstates of chaotic systems is even more 
interesting. Though the chord and centre description of the Wigner function still applies 
to mixtures of eigenstates over narrow energy windows, no classical theory accounts for 
individual eigenstates at present. Therefore, it will be extremely useful to describe 
these as interfering superpositions of gaussians placed near the energy shell. Of course, 
the size of this basis would diverge at the classical limit as — 0, but manageable 
approximations should be attainable for finite excitations. 
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